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Abstract 

In this paper, we propose a simple hybrid WENO scheme to increase computational 
efficiency and decrease numerical dissipation. Based on the characteristic-wise ap- 
proach, the scheme switches the numerical flux of each characteristic variables be- 
tween that of WENO scheme and its optimal linear scheme according to a discon- 
tinuity detector measuring the non-resolvability of the linear scheme. A number of 
numerical examples computed with 5th-order WENO scheme suggested that, while 
achieving very small numerical dissipation and good robustness, the computational 
effort on WENO flux used in the hybrid scheme is negligible. 
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1 Introduction 



High-order weighted essentially non-oscillatory (WENO) schemes 4| are non- 
linear schemes which are generally suitable for the simulation of shock-turbulence 
interaction due to their high-resolution properties. However, there are still 
two important issues hindering its extensive application. One is that, WENO 
schemes are based on local characteristic decompositions and flux splitting to 
avoid spurious oscillation. But the cost of computing nonlinear weights and 
local characteristic decompositions is very high. The other is that, WENO 
schemes are more dissipative than many low-dissipative linear schemes, par- 
ticularly in regions without strong shock wave but with large density variation 
or shear rates. 



A promising approach to overcome these drawbacks is using hybrid meth- 
ods , |8j]. The key idea is switching or blending between the non-linear 
WENO scheme and the linear scheme according to a shock sensor or disconti- 
nuity detector. Generally, there are two approaches to apply the discontinuity 

a EL 



by which the 



detector. One is the trouble-cell detecting approach 
trouble cells crossed by discontinuity are first detected by measuring the flow 
properties such as the density, pressure, entropy, etc., then the numerical fluxes 
at cell-faces of all the neighboring cells are computed by WENO scheme. Since 
the performance of this approach is strongly problem dependent, the choice 
of an effective discontinuity detector remains a problematic issue for com- 
plex applications. Another approach is applying the discontinuity detector in 



characteristic- wise 
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15] , by which the numerical fluxes of characteristic 
variables obtained by WENO scheme is switched on or blended with that of the 
linear scheme, according to the discontinuity detector usually defined with the 
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quantities corresponding to low- and high-order derivatives. Although discon- 
tinuity detector in the characteristic-wise approach is less problem dependent, 
very often multiple tunable parameters are still required to obtain more ac- 
curate or more robust solution for different problems. In addition, up to now, 
the computational efficiency in the characteristic-wise approach is either not 
considered, or is much less than that of the trouble-cell detecting approach. 

In this paper, we propose a simple high-efficiency and low-dissipation hy- 
brid WENO scheme based the characteristic-wise approach. The new scheme 
switches the numerical flux of each characteristic variables between that of 
WENO scheme and its optimal linear scheme according to a discontinuity 
detector measuring the non-resolvability of the linear scheme. Since the lin- 
ear scheme omits the computation of the WENO weights, and most of the 
characteristic-projection operations, the scheme increase the overall computa- 
tional efficiency greatly. Furthermore, by choosing the single, broadly effective 
parameter of the discontinuity detector, much lower numerical dissipation is 
achieved without compensation of robustness. 



2 Method 



We assume that the fluid is inviscid and compressible, described by the one- 
dimensional Euler equations as 

where U = (p,m,E) T , and F(U) = [m,pu 2 +p, (E + p)u} T . This set of equa- 
tions describes the conservation laws for mass density p, momentum density 
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m = pu and total energy density E = pe + pu 2 /2, where e is the internal 
energy per unit mass. To close this set of equations, the ideal-gas equation of 
state p = (7 — l)pe with a constant 7 is used. 



2.1 Characteristic- wise WENO scheme 

For completeness, we briefly recall the classical 5th-order WENO scheme A\ 
for solving Eq. ([T]). The discretization is within the spatial domain such that 
X{ = iAx, i = 0, N, where Ax is the spatial step, the semi-discretized form 
by the method of lines yields a system of ordinary differential equations 



dlL 1 



dt Ax 



■7- Fi-1/2 - Fi+1/2 , (2) 



where Fj±i/2 are the numerical flux at Xj±i/2, respectively. Once the right- 
hand side of this expression has been evaluated, TVD Runge-Kutta methods 
are employed to advance the solution in time. In the typical characteristic- 
wise finite-difference WENO scheme, the Fj + i/2 are usually approximated in 
the local characteristic fields. Let us take the matrix A i+1 / 2 to be Roe- average 
Jacobian dF/dXJ at x i+ i/ 2 . We denote by A s , r s (column vector) and \ s (row 
vector) the sth eigenvalue, right and left eigenvectors of A i+ i/ 2 , respectively. 
Then, the physical fluxes and conservative variables in the neighborhood of i 
are projected to the characteristic field by 

Uj i8 = l a -Uj, g j:S = l s ■ Fj, (3) 
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where i + 3>j>i — 2 and s = 1, 2, 3. For each component of characteristic 
variables, the corresponding split numerical fluxes are constructed by 

1 _i 

ff,s = 2 (&> + a » u i.«) ' 4- = 2 ^ s ~ asVj ^ ' 

where a s = \X S \ for a Roe flux (RF); a s = max|A/ jS |, where / represents 
the entire computational domain, for a Lax-Friedrichs flux (LF); and a s = 
max(|Aj )S |, |Aj-fi )S |), plus several other neighbors, for a local Lax-Friedrichs flux 
(LLF). The WENO reconstruction gives 

2 2 

4+1/2,8 = U k,sfk,i+l/2' fi+l/2,s = U k,sfk,i+l/2,si (5) 
fc=0 k=0 

where 4Ti±i/2 s are 3rd-order reconstruction from upwind 3-point stencils, and 
uj k s are WENO weights defined by 



U k,s - ^2 — » a M - 77, — Tq, W 



where <i = ^, rfi = | and g?2 = jq are optimal weights. These optimal 
weights generate the 5th-order upwind scheme, by which the numerical flux 
is reconstructed from a 5-point stencil, e > prevents division by zero, q — 1 
or 2 is chosen to adjust the distinct weights and k ,s are the smoothness 



indicators. We refer Jiang and Shu 



4[ for the details of WENO reconstruction. 



The numerical flux in each characteristic field is then computed by 

fi+l/2,8 — 4f.l/2,a + 4+1/2,8' (^) 

At last, this numerical flux is projected back to the component space by 

3 

Fj+l/2 = 4+l/2,s r s- (8) 
s=l 
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It can be found that most of floating-point operations in the characteristic- 
wise WENO scheme are the characteristic projection (vector-vector product) 
of Eq. ([3]), and the computation of the WENO weights in Eq. (jBJ). 



2.2 Characteristic- wise linear scheme 



In order to decrease the number of floating-point operations, we consider a 
characteristic-wise linear scheme. In this scheme, the approximated physical 
fluxes and the differences of approximated conservative variables at Xj+1/2 are 
computed component-by-component as 



Pi+1/2 = g ( F m/2 + ^+1/2) > Au m/2 = - (pt+1/2 ~ ^i+1/2) ■ ( 9 ) 



With the optimal linear scheme of the 5th-order WENO scheme, one has 



Fi+i/2 = 7^ (Fi_ 2 - 8F«_i + 37F, + 37F i+1 - 8F i+2 + F i+3 ) , (10) 
60 

AU l+1/2 = ±- (Ui_ 2 - 5Ui_! + 10U, - 10U i+1 + 5U i+2 - U l+3 ) . (11) 
60 



Note that Taylor-series expansion of Eq. (JTTj) suggests 
Ax 5 <9 5 U(x 



AU m 



/2 



60 dx 5 



+ 0(Ax 7 ), (12) 
i+1/2 



which corresponds to the 5th-order derivatives of the conservative variable. 
Then Fj+1/2 and AUj+i/2 are projected to the characteristic field by 

A%l/2, s = L ■ AU i+ i/2, 9i+l/2,s — h ■ Fj+i/2- (13) 
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For each component of the characteristic variables, the corresponding numer- 
ical flux is constructed by 



where a s are chosen in the same way as Eq. (J3J). At last, the numerical flux 
obtained in each characteristic field is projected back to the component space 
by Eq. (jSJ). Note that, compared to the characteristic- wise WENO scheme, the 
linear scheme omits the computation of the WENO weights in Eq. (EJ), and 
decreases 5/6 of characteristic-projection operations of Eq. ([3]). 

2.3 Hybrid WENO scheme 

Since the linear scheme is numerically unstable for solution with discontinu- 
ities, it should be switched on only in smooth regions of the solution. Here, by 
noticing that Av i+1 / 2)S have the dimension of density, we measure the smooth- 
ness of solution with a dimensionless discontinuity detector 



where p is the Roe- average density of A i+1 / 2 . With a given tolerance e< 1, 
the hybrid scheme is simply defined as follow: for each component of the 
characteristic variables, if a s < e, the numerical flux is obtained by linear flux 
of Eq. dH]); otherwise it is obtained by WENO flux of Eq. ©. 

It can be found from Eq. ffl2|) that cr s also correspond to the 5th-order deriva- 
tives of the characteristic variables due to the linearized characteristic projec- 
tion. Since the optimal 5th-order linear scheme reconstructs with 4th-degree 



fi+l/2,s — 9i+l/2,s + a sAv i+ i/ 2y 



(14) 




(15) 
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polynomials, and is not able to resolve 5th or higher order derivatives, a s ac- 
tually measure the non-resolvability of the linear scheme. Note that, unlike 



some discontinuity detectors 



11, 0, Ei i 



in which the quantities corresponding 
to low-order derivatives are used as the denominators, a s does not degenerate 
at critical points. Also note that the above approach can be directly applied 
to mult i- dimensions by a dimension-by-dimension approach, and can be very 
easily extended for higher-order or modified WENO schemes. The only mod- 
ification for other WENO schemes is that Eq. (j9]) should be computed with 
the corresponding optimal linear schemes. 



3 Numerical examples 



The following numerical examples are provided to illustrate the potential of the 
proposed hybrid WENO scheme. The governing equations are one- and two- 
dimensional compressible Euler or Navier-Stokes equations. While the original 
5th-order WENO is denoted as WENO-5, the present hybrid scheme is denoted 
as H- WENO-5. The 3rd-order TVD Runge-Kutta scheme is used for time 



integration [12||. The local computational efficiency on each cell-face through 
the entire computation time is measured by r\ = 1 — M WENO / 'M Linear , where 
Mweno is the number of operations with WENO flux and M Linear is the 
number of operations with linear flux. The overall computational efficiency is 
obtained by average rj over the entire computational domain. We set e = 10~ 5 
for all numerical examples (see also a study of the sensitivity of e in Sections 
13.31 and |3T3|) . If not mentioned otherwise, all the computations are carried out 
with a CFL number of 0.6. 
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3. 1 Propagation of broadband sound waves 



This problem, taken from Sun et al. 



14j . corresponds to the propagation of 



a sound wave packet which contains acoustic turbulent structure with various 
length scales. The initial condition is 



N/2 



p(x, 0) = po { 1 + e [E P (k)} 1/2 sin [2vrfc(x + <j> k )\ } , 

k=l 

/ M / (J 

p(x, 0) = p 



p(x,0) 



u(x, 0) = Mo + 



Po 
2 



c(x,0) 



Co 



where 0& is a random number between and 1 with uniform distribution, 
e = 10~ 3 , 7 = 1.4, c is the speed of sound and 



E p (k) = ( A ) exp- 2 (^) 2 



is the energy spectrum which reaches its maximum at k = kg. A periodic 
boundary condition is applied at x = and x — 1. Computations have been 
carried out on a 128-point grid using a CFL number of 0.2 for one period of 
time. The numerical results shown in Fig. [1] suggest that no WENO flux is 
switched on during the computation and the solution recovers that obtained 
by the optimal linear scheme. Note that the present methods gives consid- 
erably better resolved sound waves than WENO-5, which is numerically too 
dissipative. 
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3.2 Shock-tube problems 

Here, we show that the proposed scheme H-WENO-5 passes the shock-tube 
test problems: the Sod problem (Sod 1978), the Lax problem (Lax 1954) and 
the 123 problem (Einfeldt et al. 1991). For the Sod problem, the initial con- 
dition is 



and the final time is t — 0.2. For the Lax problem, the initial condition is 



and the final time is t — 0.14. For the 123 problem, the initial condition is 



and the final time is t = 0.1. We examine the numerical solution with three 
resolutions on grids of 200, 400 and 800 points. Figure [2] gives the density 
distributions computed on a 200-point grid, and local computational effi- 
ciency for three grid resolutions. It is observed good agreement (except with 
slightly sharp density discontinuities) between the present solution and that 
of WENO-5. Note that, as shown in Fig. [2b, d and f, the computational effi- 
ciency converges linearly to full efficiency. For all three problems, the overall 
computational efficiency achieves 99% for the lowest resolution and 99.8% for 
the highest resolution. 




(0.125,0,0.1) ifl>x>0.5 



if < x < 0.5 
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3.3 Interacting blast waves 



We consider a two-blast-wave interaction problem, which is taken from Wood- 



ward and Colella 



2] . The initial condition is 



(p,u,P) 



'(1,0,1000) if0<x<0.1 

1,0,0.01) if 0.1 < x < 0.9 
[(1,0,100) if 1 > x > 0.9 



and the final time is t — 0.038. The reflective boundary condition is applied 
at both x = and x = 1. We examine the numerical solution on a 400- 
point grid, and the reference "exact" solution is a high- resolution solution on 
a 3200-point grid computed by WENO-5. Figure [3] give the profiles of density 
and local computational efficiency Again, good agreement with the reference 
solutions is observed. Due to less numerical dissipation, H- WENO-5 shows an 
improved solution than WENO-5, especially near the valley at about x = 0.75 
and the right peak at about x = 0.78. As shown in Fig. [3b, the computational 
efficiency is only weakly dependent on the choice of e. Actually, by increasing 
e 100 times, the extra gain of overall efficiency is less than 0.5% (from 98.8% 
to 99.3%). Note that, slightly less numerical dissipation can be achieved by 
increase e (see the density profile obtained with e = 10 -4 in Fig. |3k), which 
also suggests only weak dependence. Also note that, since the computation is 
still stable after e is increased 100 times to 10 -3 , the numerical robustness of 
the present method is only weakly dependent on the choice of e too. 
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3.4 Shock- density wave interaction 



12j . The initial condi- 



We consider a shock density-wave interaction problem 
tion is set by a Mach 3 shock interacting with a perturbed density field 

(3.857,2.629,10.333) if < x < 1 
(l + 0.2sin(5x),0, 1) if 10 > x > 1 
and the final time is t — 1.8. A zero-gradient boundary condition is applied 
at x = and 10. We examine the numerical solution on a 200-point grid, 
and the reference "exact" solution is a high-resolution solution on a 3200- 
point grid computed by WENO-5. Fig. @Ji show the calculated density and 
velocity profile. It can be observed that H- WENO-5 gives considerably better 
resolved density waves behind the shock than WENO-5, which is numerically 
too dissipative. Again, as shown in Fig. lib, the computational efficiency is only 
weakly dependent on the choice of e (the overall efficiency changes from 99.3% 
to 99.6%). FiguresHb and d give the density and velocity profiles obtained with 
two different e = 5 x 10~ 5 and 10~ 4 . Again, while the numerical dissipation 
can be decreased slightly with increasing e, it is only weakly dependent on the 
later. 



3.5 Double Mach reflection a strong shock 



We consider the problem from Woodward and Colella [2| on the double Mach 
reflection of a strong shock. A Mach 10 shock in air is reflected from the wall 
with incident angle of 60°. For this case, the initial condition is 

(1.4,0,0,1) if y < 1.732(2; - 0.1667) 

(8,7.145,-4.125,116.8333) else 



(p,u,v,p) 
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and the final time is t = 0.2. The computational domain of this problem is 
[0,0] x [4,1]. Initially, the shock extends from the point x = 0.1667 at the 
bottom to the top of the computational domain. Along the bottom boundary, 
at y — 0, the region from x = to x = 0.1667 is always assigned post-shock 
conditions, whereas a reflecting wall condition is set from x = 0.1667 to x — 4. 
Inflow and outflow boundary conditions are applied at the left and right ends 
of the domain, respectively. The values at the top boundary are set to describe 
the exact motion of a Mach 10 shock. We examine the numerical solution with 
three resolutions on grids of 240 x 60, 480 x 120 and 960 x 240 points. 

Figure [5] shows contours of density and local computational efficiency com- 
puted on a 960 x 240 grid. It is observed a good agreement with the results 
of Kim and Kwon {(]] (their Figs. 12 and 13) computed with fine-tuned p_a 



1$ 



rameters at the same resolution. Compared to the results in Zhou et al 
computed with fine-tuned parameters and higher grid resolution, much more 
fine-scale structures are obtained by the present scheme. Note that, as shown 
in Fig. \5jp, WENO flux is seldom calculated in the entire computation domain, 
except near the end of the main reflection wave. The overall computational ef- 
ficiency for three resolutions are 98.3%, 99.2% and 99.6%, which again suggest 
linear convergence to full efficiency 



3.6 Viscous shock tube problem 



We consider the two-dimensional viscous flow problem in a square shock tube 
with a unit side length and insulated walls 3(. In this problem, the propagation 
of the incident shock wave and contact discontinuity lead to a thin boundary 
layer. After its reflection on the right wall, the shock wave interacts with this 
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boundary layer and results a separation region and the formation of a typical 
" A-shape like shock pattern" . The initial condition is 

;i20, 0,120/7) if < x < § 
[1.2,0,1.2/7) if 1 > ar > I 

The fluid is assumed as ideal gas with 7 = 1.4 and constant dynamics vis- 
cosity = 0.005 and Prandtl number Pr = 0.73 and satisfying the Stokes 
assumption. If the reference values are chosen as the initial speed of sound, 
unit density and unit length, the Reynolds number is 200. By applying the 
symmetry condition at the upper boundary, only the lower half domain is 
actually computed. For other boundaries, the no-slip and adiabatic wall con- 
ditions are applied. The viscous and heat transfer is calculated by 6th-order 
accuracy, Here, we exam the numerical solution with two resolutions on grids 
of 250x125 and 500x250 points. 

As shown by the density contours on the coarse grid in Fig. Eh, the computed 
result, especially on the shape and the tilting angle of the primary vortex, 
is already considerably better than those obtained by the WENO-5, the 6th- 
order artificial compression method (ACM) with wavelet based filter scheme 
131 ] (their Fig. 4) and the 5th-order mult i- dimensional limiting process (MLP) 
cocombined with M-AUSMPW+ scheme [7J (their Fig. 41) on the same grid 
resolution. Note that, numerical convergence is indicated for the results com- 
puted on the fine grid because no notable differences, except the near shock 
region, can be identified from and the converged density contours in Sjogreen 
and Yee 13J] (their Fig. 2) obtained with much higher grid resolutions. Also 
note that, the obtained overall computational efficiency for two resolutions are 
99.8%, 99.9%, which again suggested WENO flux is seldom calculated in the 
entire computation domain. 
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4 Concluding remarks 

We propose a simple hybrid WENO scheme to increase computational effi- 
ciency and decrease numerical dissipation. Based on the characteristic-wise 
approach, the scheme switches the numerical flux of each characteristic vari- 
ables between that of WENO scheme and its optimal linear scheme accord- 
ing to a discontinuity detector measuring the non-resolvability of the linear 
scheme. As shown by a number numerical examples that the overall efficiency 
are always higher than 98%, the hybrid scheme increase the computational 
efficiency greatly. Also, by choosing the single, broadly effective parameter 
of the discontinuity detector, much lower numerical dissipation is achieved 
without compensation of robustness. In addition, since the method uses a 
general characteristic-wise formulation of WENO scheme, it can be applied to 
multiple-species flows, in which the overhead of local characteristic projection 
and WENO weights calculation is even much more serious. 
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Figure 1. Propagation of broadband sound waves computed on a 128 points grid: 
(a) pressure and (b) local computational efficiency distribution. 
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Figure 4. Shock-density-wave interaction: (a) density and velocity profiles on a 
200-point grid; (b) local computational efficiency; (c, d) results obtained with 
e = 5 x 1CT 5 and 10~ 4 , respectively. 
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Figure 5. Double-Mach reflection of a Mach 10 shock wave at t = 0.2 computed on 
960 x 240 grid: (a) 40 density contours from 1.88783 to 20.9144, (b) 20 contours of 
local computational efficiency from 0.3 to 1. 
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Figure 6. Viscous shock tube problem at t = 1: 20 density contours from 15 to 125 
computed on (a) a 250 x 125 grid and (b) a 520 x 250 grid. 
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